Data mining based approach for online calibration of phasor measurement unit (pmu)

ABSTRACT

Data quality of Phasor Measurement Unit (PMU) is receiving increasing attention as it has been identified as one of the limiting factors that affect many wide-area measurement system (WAMS) based applications. In general, existing PMU calibration methods include offline testing and model based approaches. However, in practice, the effectiveness of both is limited due to the very strong assumptions employed. This invention presents a novel framework for online error detection and calibration of PMU measurement using density-based spatial clustering of applications with noise (DBSCAN) based on much relaxed assumptions. With a new problem formulation, the proposed data mining based methodology is applicable across a wide spectrum of practical conditions and one side-product of it is more accurate transmission line parameters for the energy management system (EMS) database and protective relay settings. Case studies are presented to demonstrate the effectiveness of the proposed method.

REFERENCE TO RELATED APPLICATION

This application claims the benefit of the filing date of U.S.Provisional Application Ser. No. 62,429,898, entitled “A data miningbased approach for online calibration of Phasor Measurement Unit (PMU)”and filed on Dec. 5, 2016. The teachings of the entire referencedapplication are incorporated herein by reference.

TECHNICAL FIELD

This invention is related to methods to calibrate phasor measurementunit and its applications. More particularly, this invention presents anovel framework for online error detection and calibration of PMUmeasurement using density-based spatial clustering of applications withnoise (DBSCAN) based on much relaxed assumptions.

BACKGROUND

A phasor measurement unit (PMU) is a device which measures the 50/60 HzAC waveform (voltages and currents) on an electricity grid using acommon time source Global Positioning System (GPS) radio clock forsynchronization. Time synchronization allows synchronized real-timemeasurements of multiple remote measurement points on the grid. Theresulting measurement is known as a synchrophasor. The obtainedsynchrophasor measurement is of great importance in tremendous modernpower system applications, for instance, support demand responsemechanism to manage a power system, detect early fault allowing forisolation of operative system and prevent power outages, providetrustful information for power system state estimation, etc. Phasormeasurement unit (PMU) is envisioned to be one of the enablingtechnologies in smart grid, with the promise of massive installation inthe future power systems. On one hand, most synchrophasor-basedapplications, especially the mission critical ones, require themeasurements to be very reliable and accurate. On the other hand,although PMU data are expected to be highly accurate, this potentialaccuracy and reliability are not always achieved in actual fieldinstallation due to various causes [1]. It has been observed under manyoccasions that PMU measurements can have various types of data qualityissues. To ensure accurate, reliable and consistent PMU data, there arepressing needs to calibrate PMU to fulfill the claimed performance.

As discussed in [2], the PMU device itself is typically very accurate,but the instrumental channel, where PMU gets its inputs, is usually muchless accurate. Specifically, the instrumentation channel (e.g.,potential and current transformers) can introduce magnitude and phaseangle errors that can be magnitudes of orders higher than the typicalPMU accuracy. A practically useful calibration method should be capableof handling inaccuracies originated from both PMU and itsinstrumentation channel.

Previously the Performance and Standards Task Team (PSTT) published aPMU system testing and calibration guide [3]. As discussed and widelyaccepted in the 2016 NASPI Work Group meeting, PMU data quality effortsneed to be implemented to ensure the highest synchrophasor signalquality for applications. The modified IEEE C37.118 standard requiresthe total vector error (TVE) between a measured phasor and its truevalue to be well within 1% under steady-state operating conditions [4].Towards these requirements, many PMU calibration schemes have beenproposed. In general, these methods can be divided into two categoriesbased on how they are implemented: offline testing/calibration [5-11]and model-based approaches [12-16].

The former works by comparing PMU output against standard testingsignal(s), using certain types of specialized equipment or systems whoseaccuracies are at least one level greater than the to-be-tested PMUs.This type of methods requires specialized expensive equipment/system,and due to their offline nature, errors originated from theinstrumentation channel cannot be duplicated and compensated.

The latter works by fitting PMU measurements into a mathematical modelfor fidelity check, assuming parameters of the system/device(s) and themodel are known a priori and accurate. In [12], the authors present aphasor-data-based state estimator (PSE) that is capable of identifyingand correcting bias error(s) in phase angles. This approach assumes thephasor magnitudes and network parameters are both accurate. One paper in[13] proposes the idea of a “super calibrator” for substation-level datafiltering and state estimation, the input of which includes PMU data,SCADA data, and a detailed 3-phase model of the substation, etc. Despitecomplexity of the model, the accuracy level of SCADA data addsuncertainty, or even degrades performance of the approach. Paper [1]proposes a calibration-factor-based iterative non-linear solutionapproach for 3-phase PMU data calibration. Performance of the approachis highly dependent upon accuracy of the 3-phase transmission lineparameters in the EMS database. The PMU data calibration approach in[14] again assumes the transmission line (TL) impedances are known to beexact. Papers [15] and [16] attempt to accomplish line parameterestimation and PMU calibration simultaneously, with the assumption thatone of the two PMUs generates perfect measurements, which, in practice,is really difficult to tell. The strong assumptions used in existingmodel-based methods undermine their practicability.

This invention presents a novel data mining based synchrophasormeasurement calibration framework which detects and corrects the overallsystematic or bias error(s) introduced by both PMU and itsinstrumentation channel. Major contribution of the proposed method liesin that it does not require accurate prior knowledge of the systemmathematical model/parameters. Furthermore, one byproduct of theproposed method is more accurate impedance parameters of thetransmission line for EMS database and protective relay settings. Byrelaxing those strong assumptions employed in existing model basedapproaches, the proposed method advances the practicability of onlinePMU calibration.

The remainder of this patent application is organized as follows.Section I describes the problem and related mathematical models. SectionI presents the proposed framework. Case studies are presented in sectionIII while conclusion and advantages are discussed in section IV.

SUMMARY OF THE INVENTION

This invention presents a novel approach for online calibration of PMUby using density-based spatial clustering. As compared to existingmethods, this invention has two major advantages: 1) it identifies theoverall bias errors introduced by both PMU and its instrumentationchannel; 2) it does not require accurate system model/parameters.Therefore, it is applicable across a wide spectrum of practicalconditions. In addition, one by-product of the proposed approach is moreaccurate transmission line (TL) impedance estimation for improved systemmodeling, more accurate protective relay settings, and other relatedapplications.

BRIEF DESCRIPTIONS OF FIGURES

FIG. 1 Bias error in PMU measurement.

FIG. 2 Transmission line nominal/equivalent PI model.

FIG. 3 Sensitivity analysis-estimated bias error in phase angle vs. linereactance.

FIG. 4 Sensitivity analysis-estimated bias error in magnitude vs. linereactance.

FIG. 5 Schematic diagram of DB SCAN with minPts=4.

FIG. 6 Feasible region for transmission line (TL) impedances.

FIG. 7 An example of bias error identification using DB SCAN.

FIG. 8 Flow chart of proposed PMU data calibration approach.

FIG. 9 Relationships between identified bias errors and error inR_(EMS).

FIG. 10 Results of DBSCAN for Case II.

FIG. 11 DBSCAN Clustering Results for Case III.

FIG. 12 Experiment result with real PMU data (errors in R, X, and Bc vs.cluster size).

DETAILED DESCRIPTION OF THE ILLUSTRATED EMBODIMENTS I. ProblemDescription and Formulation A. Bias Error in PMU Measurements

Generally speaking, errors in synchrophasor measurements can originatefrom three possible sources as discussed in [17]: transducers,synchronization, and phasor estimation algorithm. Impacts of these threesources can be summarized into two categories: random error and bias(systematic) error.

Random error, as its name suggests, is random in either direction in itsnature and difficult to predict. Random error can be circumvented frommeasurements via statistical means. Extensive studies have beenconducted in reducing random error or its influences to PMUmeasurements, with satisfactory results observed: unbiased linear leastsquares (LS) is used in [18]; non-linear LS algorithms are used in [19,20]; total LS is introduced in [21]; other optimization procedures arediscussed in [22, 23].

TABLE I MAXIMUM BIAS ERRORS INTRODUCED BY DIFFERENT PORTIONS OF THE PMUMEASUREMENT CHAIN IN A 230-KV SYSTEM [2] Voltage Phasor Current PhasorMagnitude Phasor Angle Magnitude Phasor Angle Source (%) (degree) (%)(degree) of PT/CT Error 0.6 1.04 1.2 0.52 Cabling 100 ft 400 ft 100 ft400 ft 100 ft 400 ft 100 ft 400 ft 0.009 0.009 0.115 0.411 0.066 0.0660.03 0.09 PMU 0.1 0.02 0.1 0.02 Total 0.709 0.709 1.175 1.471 1.3661.366 0.57 0.63

Systematic or bias error is reproducible inaccuracy that is consistentin the same direction. Bias error is much harder to estimate and remove.Authors of [2] have examined the maximum bias errors introduced bydifferent portions of the measurement chain. Table I summarizes themaximum bias error for a typical 230-kV system. For example, with a400-ft instrumentation cable, the maximum bias errors in the magnitudeand phase angle of the voltage phasor are 0.709% and 1.471 degree,respectively. These bias errors are no longer negligible and asystematic approach needs to be developed to identify and remove them,which is the scope of this invention.

B. Notations and Models

FIG. 1 shows a measured voltage phasor V, its corresponding true phasorvalue V _(true), and the associated bias errors in magnitude ∂V andphase angle ∂θ. The following relationship is derived:

V=V·e ^(jθ) ^(V)   (1)

V _(true)=(V+∂V)·e^(j(θ) ^(V) ^(+∂θ))  (2)

where V and θ_(V) are the magnitude and phase angle of phasor V,respectively.

A PI model as shown in FIG. 2 is considered in this work as the modelfor a general transmission line. It could either be a nominal PI modelif the line is short or an equivalent if the line is longer [24]. inFIG. 2, V _(s) and Ī_(s) represent the positive sequence voltage andcurrent phasors measured at sending end of the line while V _(r) andĪ_(r) are the corresponding phasors collected from receiving end.Variables Z and Y represent series impedance and shunt admittance of thePI model.

The following equations are derived from nodal analysis:

$\begin{matrix}{{{\overset{\_}{I}}_{s} - {{\overset{\_}{V}}_{s} \cdot \frac{Y}{2}} + {\overset{\_}{I}}_{r} - {{\overset{\_}{V}}_{r} \cdot \frac{Y}{2}}} = 0} & (3) \\{{{\overset{\_}{V}}_{s} - {Z \cdot \left( {{\overset{\_}{I}}_{s} - {{\overset{\_}{V}}_{s} \cdot \frac{Y}{2}}} \right)} - {\overset{\_}{V}}_{r}} = 0} & (4)\end{matrix}$

and

Z=R+jX  (5)

Y=G+jB _(c)  (6)

where G and B_(c) are line shunt conductance and susceptance.Combing equations (3)-(4) to solve for Y and Z yields:

$\begin{matrix}{Z = \frac{{\overset{\_}{V}}_{s}^{2} - {\overset{\_}{V}}_{r}^{2}}{{{\overset{\_}{I}}_{s} \cdot {\overset{\_}{V}}_{r}} - {{\overset{\_}{I}}_{r} \cdot {\overset{\_}{V}}_{s}}}} & (7) \\{Y = {2 \cdot \frac{{\overset{\_}{I}}_{s} + {\overset{\_}{I}}_{r}}{{\overset{\_}{V}}_{s} + {\overset{\_}{V}}_{r}}}} & (8)\end{matrix}$

Substituting each phasor in (7)-(8) with its magnitude and phase angleaccording to (1) and setting phase angle of I_(r), θ_(I) _(r) , as thereference, the following equations can be derived:

$\begin{matrix}{R = {{real}\left( \frac{{V_{s}^{2}e^{i\; 2\; \theta_{V_{s}}^{\prime}}} - {V_{r}^{2}e^{i\; 2\; \theta_{V_{r}}^{\prime}}}}{{{I_{s} \cdot V_{r}}e^{i({\theta_{I_{s}}^{\prime} + \theta_{V_{r}}^{\prime}})}} - {{I_{r} \cdot V_{s}}e^{i\; \theta_{V_{s}}^{\prime}}}} \right)}} & (9) \\{X = {{imag}\left( \frac{{V_{s}^{2}e^{i\; 2\; \theta_{V_{s}}^{\prime}}} - {V_{r}^{2}e^{i\; 2\; \theta_{V_{r}}^{\prime}}}}{{{I_{s} \cdot V_{r}}e^{i({\theta_{I_{s}}^{\prime} + \theta_{V_{r}}^{\prime}})}} - {{I_{r} \cdot V_{s}}e^{i\; \theta_{V_{s}}^{\prime}}}} \right)}} & (10) \\{B_{c} = {{imag}\left( {2 \cdot \frac{{I_{s}e^{i\; \theta_{I_{s}}^{\prime}}} + I_{r}}{{V_{s}e^{i\; \theta_{V_{s}}^{\prime}}} + {V_{r}e^{i\; \theta_{V_{r}}^{\prime}}}}} \right)}} & (11)\end{matrix}$

where θ′_(V) _(s) =θ_(V) _(s) −θ_(I) _(r) , θ′_(V) _(r) =θ_(V) _(r)−θ_(I) _(r) , and θ′_(I) _(s) =θ_(I) _(s) −θ_(I) _(r) .Shunt conductance of a transmission line, G, is usually very small andtherefore neglected from the PI model.

C. Least Square Estimator (LSE)

To investigate the sensitivity of line parameters to bias errors in PMUmeasurements, partial derivatives need to be taken for (9)-(11), all ofwhich are complex equations. For the derivatives to be valid, they mustobey the Cauchy-Riemann equations [25]. The compliancechecking/procedure is not discussed here due to space consideration, butthe validation has been completed. The following equations have beenderived:

∂R=A _(R) ·∂V _(s) +B _(R) ·∂V _(r) +C _(R) ·∂I _(s) +D _(R) ·∂I _(r) +E_(R)··∂θ′_(V) _(s) +F _(R)·∂θ′_(V) _(r) +G _(R)·∂θ′_(I) _(s)   (12)

∂X=A _(X) ·∂V _(s) +B _(X)·∂V _(r) +C _(X)·∂I_(s) +D _(X) ·∂I _(r) +E_(X)·∂θ′_(V) _(s) +F _(X)·∂θ′_(V) _(r) +G _(X)·∂θ′_(I) _(s)   (13)

∂B _(c) =A _(B) ·∂V _(s) +B _(B) ·∂V _(r) +C _(B) ·∂I _(s) +D _(B) ·∂I_(r) +E _(B)·∂θ′_(V) _(s) +F _(X)·∂θ′_(V) _(r) +G _(X)·∂θ′_(I) _(s)  (14)

where coefficients A_(x)˜G_(x) are all partial derivatives. Taking R asan example, these coefficients are:

${A_{R} = \frac{\partial R}{\partial V_{s}}},{B_{R} = \frac{\partial R}{\partial V_{r}}},{C_{R} = \frac{\partial R}{\partial I_{s}}},{D_{R} = \frac{\partial R}{\partial I_{r}}},{E_{R} = \frac{\partial R}{\partial\theta_{V_{s}}^{\prime}}},{F_{R} = \frac{\partial R}{\partial\theta_{V_{r}}^{\prime}}},{{{and}\mspace{14mu} G_{R}} = {\frac{\partial R}{\partial\theta_{I_{s}}^{\prime}}.}}$

For space consideration, detailed information of these partialderivatives is discussed in Appendix I. Put equations (12)-(14) intomatrix form to obtain:

$\begin{matrix}{\begin{bmatrix}{\partial R} \\{\partial X} \\{\partial B_{c}}\end{bmatrix} = {\begin{bmatrix}A_{R} & B_{R} & C_{R} & D_{R} & E_{R} & F_{R} & G_{R} \\A_{X} & B_{X} & C_{X} & D_{X} & E_{X} & F_{X} & G_{X} \\A_{B} & B_{B} & C_{B} & D_{B} & E_{B} & F_{B} & G_{B}\end{bmatrix} \cdot \begin{bmatrix}{\partial V_{s}} \\{\partial V_{r}} \\{\partial I_{s}} \\{\partial I_{r}} \\{\partial\theta_{V_{s}}^{\prime}} \\{\partial\theta_{V_{r}}^{\prime}} \\{\partial\theta_{I_{s}}^{\prime}}\end{bmatrix}}} & (15)\end{matrix}$

It should be noted that coefficients A_(x)˜G_(x) vary with the loading(current), as can be seen from the expression of, for example,

${C_{R} + {iC}_{X}} = {\frac{\partial Z}{\partial I_{s}} = {\frac{V_{r}{e^{i({\theta_{V_{r}}^{\prime} + \theta_{I_{s}}^{\prime}})} \cdot \left( {{V_{r}^{2}e^{i\; 2\theta_{V_{r}}^{\prime}}} - {V_{s}^{2}e^{i\; 2\theta_{V_{s}}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{V_{s}}^{\prime}}} - {I_{s}V_{r}e^{i\; {({\theta_{V_{r}}^{\prime} + \theta_{I_{s}}^{\prime}})}}}} \right)^{2}}.}}$

Assuming N sets of PMU measurements under different load conditions arecollected, the following matrix can be written:

$\begin{matrix}{H = \begin{bmatrix}A_{R\; 1} & B_{R\; 1} & C_{R\; 1} & D_{R\; 1} & E_{R\; 1} & F_{R\; 1} & G_{R\; 1} \\A_{X\; 1} & B_{X\; 1} & C_{X\; 1} & D_{X\; 1} & E_{X\; 1} & F_{X\; 1} & G_{X\; 1} \\A_{B\; 1} & B_{B\; 1} & C_{B\; 1} & D_{B\; 1} & E_{B\; 1} & F_{B\; 1} & G_{B\; 1} \\A_{R\; 2} & B_{R\; 2} & C_{R\; 2} & D_{R\; 2} & E_{R\; 2} & F_{R\; 2} & G_{R\; 2} \\A_{X\; 2} & B_{X\; 2} & C_{X\; 2} & D_{X\; 2} & E_{X\; 2} & F_{X\; 2} & G_{X\; 2} \\A_{B\; 2} & B_{B\; 2} & C_{B\; 2} & D_{B\; 2} & E_{B\; 2} & F_{B\; 2} & G_{B\; 2} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\A_{RN} & B_{RN} & C_{RN} & D_{RN} & E_{RN} & F_{RN} & G_{RN} \\A_{XN} & B_{XN} & C_{XN} & D_{XN} & E_{XN} & F_{XN} & G_{XN} \\A_{BN} & B_{BN} & C_{BN} & D_{BN} & E_{BN} & F_{BN} & G_{BN}\end{bmatrix}} & (16) \\{{F = \begin{bmatrix}{\partial V_{s}} \\{\partial V_{r}} \\{\partial I_{s}} \\{\partial I_{r}} \\{\partial\theta_{V_{s}}^{\prime}} \\{\partial\theta_{V_{r}}^{\prime}} \\{\partial\theta_{I_{s}}^{\prime}}\end{bmatrix}},{E = \begin{bmatrix}{\partial R_{1}} \\{\partial X_{1}} \\{\partial B_{c\; 1}} \\{\partial R_{2}} \\{\partial X_{2}} \\{\partial B_{c\; 2}} \\\vdots \\{\partial R_{n}} \\{\partial X_{n}} \\{\partial B_{cn}}\end{bmatrix}}} & {(17)\text{-}(18)}\end{matrix}$

If an accurate set of line impedance parameters is known a priori, thebias errors in the PMU measurements can be easily estimated using thestandard least square estimator, as:

F=(H ^(T) H)⁻¹ H ^(T) E  (19)

Seven unknowns appear in F and therefore the rank of H matrix has to beno less than seven which requires 3×N≥7 or N≥3, (N∈N*). Vector E iscomprised of the difference between the true line impedance values andthe calculated ones using (9)-(11). With the assumption that accurateline impedances are known, here are the steps for evaluating the biaserrors in PMU measurements:

Step 1: calculate line impedance parameters R, X, and B according to(9)-(11);

Step 2: evaluate vector E by comparing the calculated line impedances(from step 1) to their corresponding references obtained from the EMSdatabase: R_(EMS), X_(EMS), and B_(EMS);

Step 3: evaluate matrix H with the partial derivatives calculated fromPMU measurements: V _(s), V _(r), Ī_(s), and Ī_(r);

Step 4: solve for vector F based on (19).

D. The Challenge

The aforementioned least square estimator is able to identify biaserrors in PMU measurements assuming the line's actual impedances, as thereferences, are known a priori. In practice, these parameters are readoff from EMS database and were originally calculated based on towergeometries, conductor dimensions, estimates of line length, andconductor sags, etc. They only approximate the effects of conductor sagsand ignore the dependence of impedance parameters on temperature andloading conditions [26, 27]. Therefore, the challenge is that onlyapproximates of line impedances are known and without knowing their truevalues the calculated bias errors might be far from being accurate. Inthe following section, this invention shows how this challenge can beaddressed by using data mining technology like Density-based SpatialClustering.

II. Proposed Solution A. Sensitivity Analysis

This subsection conducts sensitivity analysis and investigates influenceof errors in referenced line impedances on bias error estimation in PMUmeasurements. A simulated transmission line with specifications shown inAppendix II is used for this study. Errors are added to all three lineimpedance references, one at a time, and the least square estimatordescribed in section I. C. is employed to evaluate F. As an example,results of the sensitivity analysis for line reactance X are shown inFIG. 3-FIG. 4.

From FIG. 3 FIG. 4, it can be observed that the influence of error inX_(EMS) on bias error estimation is linear. A 10% error in X_(EMS) willresult in: 0.16 rad error in ∂θ′_(V) _(s) and ∂θ′_(V) _(r) ; 0.01 raderror in ∂θ′_(I) _(s) ; 0.66 pu error in ∂I_(r); 0.28 pu error in∂I_(s); 0.20 pu error in both ∂V_(s) and ∂V_(r). Bias error estimationin θ′_(I) _(s) is not really affected by error in X_(EMS) while itscorresponding magnitude is affected the most among all bias errors.Results of the sensitivity analysis are summarized in Table II.

TABLE II RESULTS FOR SENSITIVITY ANALYSIS Error ∂V_(s) ∂V_(r) ∂I_(s)∂I_(r) ∂θ′_(Vs) ∂θ′_(Vr) ∂θ′_(Is) (%) p.u. rad ∂R −10 0.051 0.051 0.0720.152 0.005 0.005 0.001 −5 0.025 0.025 0.036 0.075 0.003 0.003 0.000 00.000 0.000 0.000 0.000 0.000 0.000 0.000 ∂X −10 0.200 0.200 0.283 0.6580.163 0.162 0.010 −5 0.101 0.101 0.142 0.329 0.082 0.081 0.005 0 0.0000.000 0.000 0.000 0.000 0.000 0.000 ∂BC −10 0.061 0.061 0.162 0.3600.090 0.080 0.010 −5 0.030 0.030 0.081 0.181 0.045 0.040 0.005 0 0.0000.000 0.000 0.000 0.000 0.000 0.000

A few observations can be made based on Table II, 1) PMU measurementbias error(s) estimation is generally sensitive to error(s) in lineimpedance references; 2) impact of reference errors to bias errorestimation is linear; 3) if line impedance references are exactly knowna priori, all bias error can be accurately estimated.

However, in practice, line impedance references cannot be known exactly,as discussed in section I. D. To address this challenge, a data miningapproach based on DBSCAN is proposed to address the uncertainties intransmission line (TL) impedance references.

B. DBSCAN Basics

Density-based spatial clustering of applications with noise (DBSCAN) isan unsupervised data mining technique which is able to classify datapoints of any dimension into core points, reachable points and outliers[28]. A core point p contains at least minPts points (including p)within the designated searching distance ϵ. A reachable point q existsif there exists a path p1, p2, . . . , q, so that all points on thepath, except q, are core points. Points that are not reachable from anyother point are outliers. Core points and reachable points can form acluster while outliers are excluded from such cluster. FIG. 5 shows anexample for DBSCAN with minPts=4. Note that setting minPts to 3 willgenerate the same clustering result.

As shown in FIG. 5, core points are in red, each of which has at least 4points with distance less than ϵ. The yellow ones (reachable points) arereachable from the red ones but do not have the required minimum numberof points nearby within the distance of ϵ. The blue one is not reachablefrom any other point and therefore is an outlier. The red and yellowpoints form a cluster with the blue one excluded.

C. PMU Calibration Using DBSCAN

Although EMS references can be significantly wrong, our experience showsthat the error bands are generally well within 20%. Therefore, we maydefine a as the error band multiplier for impedance references obtainedfrom the EMS database (R_(EMS), X_(EMS), and B_(EMS)). The followingconstraints can be considered:

$\begin{matrix}\left\{ \begin{matrix}{{\left( {1 - \alpha} \right)R_{EMS}} \leq R \leq {\left( {1 + \alpha} \right)R_{EMS}}} \\{{\left( {1 - \alpha} \right)X_{EMS}} \leq X \leq {\left( {1 + \alpha} \right)X_{EMS}}} \\{{\left( {1 - \alpha} \right)B_{EMS}} \leq B_{c} \leq {\left( {1 + \alpha} \right)B_{EMS}}}\end{matrix} \right. & (20)\end{matrix}$

The corresponding feasibility region can be visualized as the cube shownin FIG. 6.

The basic idea of the invented approach is to 1) scan every point withinfeasible region (a total of M points); 2) evaluate the correspondingbias errors in the PMU measurements; 3) form sets of points with eachset containing seven 4-dimensional data points, and each data point hasthe form of (∂R, ∂X, ∂B_(c), x), where x is one of the bias errors inPMU measurements (a total of M sets). 4) apply DBSCAN to cluster all theM data sets to find out the one with least number of outliers (maximumnumber of core and reachable points) and minimum searching distance.Once this cluster is identified, the actual bias errors in all PMUchannels and errors in line impedance references can be determinedaccordingly.

To minimize the computation, equation (19) is extended to (21). Ascompared to vector E and F in (19), matrix E′(3N-by-M) and F′(7-by-M)are the extended version which relates multiple bias error sets tomultiple sets of the error in referenced impedances. FIG. 7 shows asimple illustrative example of the DBSCAN algorithm for bias erroridentification.

$\mspace{745mu} {{(21)\begin{bmatrix}{\partial V_{s}^{1}} & {\partial V_{s}^{2}} & \; & {\partial V_{s}^{M}} \\{\partial V_{r}^{1}} & {\partial V_{r}^{2}} & \; & {\partial V_{r}^{M}} \\{\partial I_{s}^{1}} & {\partial I_{s}^{2}} & \ldots & {\partial I_{s}^{M}} \\{\partial I_{r}^{1}} & {\partial I_{r}^{2}} & \ldots & {\partial I_{r}^{M}} \\{{\partial\theta}\; V_{s}^{\prime 1}} & {{\partial\theta}\; V_{s}^{\prime 2}} & \ldots & {{\partial\theta}\; V_{s}^{\prime \; M}} \\{{\partial\theta}\; V_{r}^{\prime 1}} & {{\partial\theta}\; V_{r}^{\prime 2}} & \; & {{\partial\theta}\; V_{r}^{\prime \; M}} \\{{\partial\theta}\; I_{s}^{\prime 1}} & {{\partial\theta}\; I_{s}^{\prime 2}} & \; & {{\partial\theta}\; I_{s}^{\prime \; M}}\end{bmatrix}} = {\left( {H^{T}H} \right)^{- 1}{H^{T}\begin{bmatrix}{\partial R_{1}^{1}} & {\partial R_{1}^{2}} & \; & {\partial R_{1}^{M}} \\{\partial X_{1}^{1}} & {\partial X_{1}^{2}} & \; & {\partial X_{1}^{M}} \\{\partial B_{c\; 1}^{1}} & {\partial B_{c\; 1}^{2}} & \; & {\partial B_{c\; 1}^{M}} \\{\partial R_{2}^{1}} & {\partial R_{2}^{2}} & \; & {\partial R_{2}^{M}} \\{\partial X_{2}^{1}} & {\partial X_{2}^{2}} & \ldots & {\partial X_{2}^{M}} \\{\partial B_{c\; 2}^{1}} & {\partial B_{c\; 2}^{2}} & \ldots & {\partial B_{c\; 2}^{M}} \\\vdots & \vdots & \ldots & \vdots \\\vdots & \vdots & \; & \vdots \\\vdots & \vdots & \; & \vdots \\{\partial R_{N}^{1}} & {\partial R_{N}^{2}} & \; & {\partial R_{N}^{M}} \\{\partial X_{N}^{1}} & {\partial X_{N}^{2}} & \; & {\partial X_{N}^{M}} \\{\partial B_{cN}^{1}} & {\partial B_{cN}^{2}} & \; & {\partial B_{cN}^{M}}\end{bmatrix}}}}$

Or

F′=(H ^(T) H)⁻¹ H ^(T) E′

FIG. 7 shows an exemplary DBSCAN result for one data set (one column inF′). Point zero is pre-defined as one of the core point and the clustersearch always starts from zero. According to IEEE standard [4], thesearching distance E can be initially set to 0.03% for phasor magnitudeand 0.01° for phasor angle. The minimum number requirement minPts is setto be 3. In this particular example, the clustering results indicatesthat a total of 6 points are either core or reachable points while 2points are outliers (the ones with d∂I_(s) and ∂θ′_(I) _(s) ). Thepseudo-code for the proposed DB SCAN based PMU calibration is presentedbelow.

algorithm DBSCAN_based_PMU_calibration is    input: PMU phasormeasurements,     Nodes (R, X, B_(c)) in the cube.    output: PMU biaserrors      matrix(H) ←PMU measurements      matrix(E′) ←Matrix(Nodes(R,X, B_(c)) - EMS(R, X, B_(c)))   //(LSE: Least Square Estimator)     matrix(F′) ←LSE(H, E′)   vector(CCP)← //CCP: count cluster point  minPts← 3    for each column in matrix(F′) do    // i← row index formatrix (F′), (1≤i≤7)    // j← column index for matrix (F′)     residual(i, j)← |F′ (i, j)| (1≤i≤ 7)      if residual(i,j) iswithin the boundary α       C(i, j) ← F′(i, j)       CCP(j)++   Endforeach    for each column in matrix(C) do      // i← row indexfor matrix (F′), (1≤i≤ 7)   // j← column index for matrix (C)     distance(i, i ± 1)← |C(i, j) − C(i ± 1, j)|. //find the distancefrom each cluster point to its neighbor cluster point      eps(j)=(maximum(|distance|),minPts) // eps(j) is the cluster searching distancefor j^(th) column  Endforeach  k← column in F′ with maximum CPP andminimum eps    output

A flowchart of the proposed data mining based PMU data calibrationapproach is shown in FIG. 8. As shown in the sensitivity analysis ofsection II, impacts of EMS reference errors to the calculation of biaserror in each PMU measurement are given. In the ideal case when there isno error in the EMS reference and no bias error in the PMU measurements,by using DBSCAN, all the seven curves in the 4-dimensional space (∂R,∂X, ∂Bc, x), where x is one of the bias errors in PMU measurements,intersect at a single point (0, 0, 0, 0). When there are errors eitherin the EMS references or in the PMU measurements, according to theproposed approach, coordinates of this intersection corresponds to theerror(s) in the EMS reference, while the outlier(s) identifiedcorrespondingly are the bias errors in PMU measurements.

The aforementioned discussion is for the situation without noise. Whennoise exists in PMU measurements, the seven curves typically will notintercept exactly at a single point but instead will stay very close toeach other around one particular zone/region, which may be referred toas the “zero region”. By looking at the searching distance and number ofoutliers, the “zero region” can be identified and therefore the errorsin EMS references and PMU measurements can be evaluated accordingly.

III. Experimental Validation

Five case studies are presented in this invention to demonstrate theprocedure and effectiveness of the proposed PMU data calibrationframework. A testing system with parameters shown in Table VIII has beenset up in Matlab/Simulink for these experiments.

A. Case I: The Ideal Case-No Error in Impedance References

Objective of the first case study is to validate performance of theproposed method when no error exists in the TL impedance references.Basically, different sets of combinations of bias errors have been addedto the PMU measurements and the proposed approach is used to identifythem. The results for six representative cases are summarized in TableIII, in which both the true bias errors and the calculated ones arepresented and compared. The agreement between true bias errors andcalculated ones validates the proposed approach under the idealcondition with no error in the referenced impedances.

TABLE III BIAS ERROR IDENTIFICATION UNDER THE IDEAL CONDITION TrueCalculated True Calculated True Calculated p.u. or rad (×10⁻³) ∂V_(s) 00.0334 10 9.9634 10 10.0395 ∂V_(r) 0 0.0335 10 9.9613 10 10.0402 ∂I_(s)10 10.0591 0 −0.2240 10 10.0177 ∂I_(r) 0 0.1480 0 −0.5332 0 −0.2318∂θV_(s)′ 0 0.0067 1.75 1.7836 0 0.0185 ∂θV_(r)′ 1.75 1.7601 0 0.06031.75 1.7536 ∂θI_(s)′ 0 −0.0001 0 −0.0070 0 −0.0070 ∂V_(s) 10 10.0161 1010.0891 10 10.0033 ∂V_(r) 10 10.0154 10 10.0894 10 10.0034 ∂I_(s) 109.9647 10 10.0591 10 9.9328 ∂I_(r) 0 −0.3783 10 9.9659 10 9.6749∂θ′_(Vs) 0 0.0579 0 0.0138 1.75 1.7731 ∂θ′_(Vr) 1.75 1.7937 1.75 1.75221.75 1.7757 ∂θ′_(Is) 1.75 1.7663 1.75 1.7698 1.75 1.7676

B. Case II: One Referenced Impedance Has Error

The second case study considers error in one of the referencedimpedances. Towards this goal, errors are added to each of thereferenced impedances, one at a time, and different combinations of biaserrors are considered for PMU measurements. For space consideration,only the result for a representative case is presented below. And inthis particular case, a −2% error is considered for the seriesresistance, R_(EMS), and 0.01 p.u. bias error is added to magnitude ofthe sending-end current phasor. A 20% error band is considered for thereferenced impedances with a being set to 20%. The proposed approachscans all 4-dimensional data points collected from matrix F′, and forvisualization purpose, only the relationship between bias errors anderrors in R is plotted as shown in FIG. 9. The dashed line marks theoutcome of DB SCAN and the X-axis gives the corresponding error inR_(EMS). To help illustrate the DBSCAN process, FIG. 10 visualizes thecorresponding DBSCAN clustering results, in terms of size of the clusterand searching distance. Final results are summarized in Table IV.

TABLE IV TEST RESULTS FOR CASE II Calculated Error in True Calculated RX B_(C) p.u. or rad (×10⁻³) (%) (%) (%) ∂V_(s) 0 0.0230 −2.02 0.01 −0.01∂V_(r) 0 0.0345 ∂I_(s) 10 10.0328 ∂I_(r) 0 −0.0175 ∂θ′_(Vs) 0 0.0224∂θ′_(Vr) 0 0.0130 ∂θ′_(Is) 0 −0.0149

According to Table IV, the proposed approach successfully identifies notonly bias error in PMU measurements but also error in the referenced TLseries resistance.

C. Case III. Two Referenced Impedances Have Errors

In the third case study, −4% error and −6% error are considered forR_(EMS) and X_(EMS), respectively; bias errors of 0.01 p.u. and 0.00175rad are added to V_(s) and θ_(Vr), respectively. A 20% error band isconsidered for the referenced impedances with α being set to 20%.

To help illustrate the DBSCAN process, FIG. 11 visualizes the clusteringresults, in terms of size of the cluster and searching distance. Botherrors in the referenced impedances and the bias errors in PMUmeasurements are successfully identified. Final results are summarizedin Table V.

TABLE V TEST RESULTS FOR CASE III Calculated Error in True Calculated RX B_(C) p.u. or rad (×10⁻³) (%) (%) (%) ∂V_(s) 10 10.0612 −4.2 −5.8 0.02∂V_(r) 0 0.0334 ∂I_(s) 0 0.0316 ∂I_(r) 0 0.1591 ∂θ′_(Vs) 0 0.0036∂θ′_(Vr) 1.75 1.7617 ∂θ′_(Is) 0 0.0047

D. Case IV. All Referenced Impedance Values Have Errors

In the fourth case study, a set of −2%, −5%, 2% errors are consideredfor R_(EMS), X_(EMS), and B_(EMS), respectively; bias errors of 0.01p.u. and 0.00175 rad are added to V_(s) and θV_(r), respectively. A 20%error band is considered for the referenced impedances with α being setto 20%. Experimental results are shown in Table VI. Table VI demonstrateagain the effectiveness of the proposed method when all referencedimpedance values have errors. One key value of the proposed approachlies in its capability of PMU calibration without knowing an accuratesystem model.

TABLE VI TEST RESULTS FOR CASE IV Calculated Error in True Calculated RX B_(C) p.u. or rad (×10⁻³) (%) (%) (%) ∂V_(s) 10 10.0143 −2.1 −5.2 2.2∂V_(r) 0 0.0054 ∂I_(s) 0 −0.0375 ∂I_(r) 0 −0.0846 ∂θ′_(Vs) 1.75 1.7967∂θ′_(Vr) 0 −0.0437 ∂θ′_(Is) 0 −0.0091

E. Case V. Experiment Using Real PMU Data

In the fifth case study, PMU data are collected for a 500-kVtransmission line in Jiangsu Electricity Power Grid with the name of“Huadong-Tianhui Line #5621”. PMU data reporting rate is 25 samples persecond.

Using these real PMU data, the proposed approach is applied to identifyboth errors from EMS database and bias errors in the measured phasors.As discussed above, cluster size and searching distance, ϵ, are used asthe criteria for DBSCAN. Part of the spatial clustering results arevisualized in FIG. 12, which shows the relationship between errors inthe reference values of R, X, and B_(c), and the cluster size. As shownin FIG. 12, the top three points with the largest cluster size areA1˜A3, of which A2 has the smallest searching distance. The finalresults are summarized in Table VII. The results show:

1) TL impedance parameters, R, X and Bc in the EMS database have errorsof −14%, 5.4% and 12.6% respectively;

2) no significant bias error in the voltage phasors collected from thereal PMU is identified;

3) bias errors of 0.0171 pu and 0.0164 pu are identified in themagnitudes of sending-end and receiving-end current phasors,respectively;

4) no significant bias error is identified in the phase angles of thetwo current phasors.

TABLE VII TEST RESULTS FOR CASE V Calculated Error in Calculated R XB_(C) p.u. or rad (×10⁻³) (%) (%) (%) ∂V_(s) −0.0032 −14.0 6.4 12.6∂V_(r) −0.0033 ∂I_(s) 0.0171 ∂I_(r) 0.0164 ∂θ′_(Vs) 0.0003 ∂θ′_(Vr)0.0027 ∂θ′_(Is) −5.6238e−5

Computation time of the proposed approach is dependent upon the numberof data points to scan within the feasibility region. In one experiment,a total of 1 million points (worst case scenario) are processed using aMatlab program, and the solution process takes roughly 29 seconds(recoding the program using C++ will greatly speed up the solutionprocess). Fortunately, PMU data calibration does not need to beconducted very often, and once a week or longer will work for mostcases. Structure of the proposed algorithm is suitable for parallelprocessing, which will further speed up the solution process.

As show in the above diagram, this invention can be implemented inhardware, firmware or software, or a combination of the three.Preferably the invention is implemented in a computer program executedon a programmable computer having a processor, a data storage system,volatile and non-volatile memory and/or storage elements, at least oneinput device and at least one output device.

By way of example, a block diagram of a computer to support the systemis discussed next. The computer preferably includes a processor, randomaccess memory (RAM), a program memory (preferably a writable read-onlymemory (ROM) such as a flash ROM) and an input/output (I/O) controllercoupled by a CPU bus. The computer may optionally include a hard drivecontroller which is coupled to a hard disk and CPU bus. Hard disk may beused for storing application programs, such as the present invention,and data. Alternatively, application programs may be stored in RAM orROM. I/O controller is coupled by means of an I/O bus to an I/Ointerface. I/O interface receives and transmits data in analog ordigital form over communication links such as a serial link, local areanetwork, wireless link, and parallel link. Optionally, a display, akeyboard and a pointing device (mouse) may also be connected to I/O bus.Alternatively, separate connections (separate buses) may be used for I/Ointerface, display, keyboard and pointing device. Programmableprocessing system may be preprogrammed or it may be programmed (andreprogrammed) by downloading a program from another source (e.g., afloppy disk, CD-ROM, or another computer).

Each computer program is tangibly stored in a machine-readable storagemedia or device (e.g., program memory or magnetic disk) readable by ageneral or special purpose programmable computer, for configuring andcontrolling operation of a computer when the storage media or device isread by the computer to perform the procedures described herein. Theinventive system may also be considered to be embodied in acomputer-readable storage medium, configured with a computer program,where the storage medium so configured causes a computer to operate in aspecific and predefined manner to perform the functions describedherein.

The invention has been described herein in considerable detail in orderto comply with the patent Statutes and to provide those skilled in theart with the information needed to apply the novel principles and toconstruct and use such specialized components as are required. However,it is to be understood that the invention can be carried out byspecifically different equipment and devices, and that variousmodifications, both as to the equipment details and operatingprocedures, can be accomplished without departing from the scope of theinvention itself.

These are key steps of this invention:

-   -   1. Least square estimator construction based on PI model, PMU        measurements and TL parameters from EMS;    -   2. Set confidential interval to EMS data;    -   3. Estimation result data clustering using DB SCAN;    -   4. Cluster filter targeting the answer indicates bias errors in        PMU measurements.

More detailed steps of this invention:

-   -   1) scan every point within feasible region (a total of M        points);    -   2) evaluate the corresponding bias errors in the PMU        measurements;    -   3) form sets of points with each set containing seven        4-dimensional data points, and each data point has the form of        (∂R, ∂X, ∂B_c, x), where x is one of the bias errors in PMU        measurements (a total of M sets);    -   4) apply DBSCAN to cluster all the M data sets to find out the        one with least number of outliers (maximum number of core and        reachable points) and minimum searching distance;    -   5) determine actual bias errors in all PMU channels and errors        in line impedance references after the cluster is identified in        step 4).

The following flow chart shows detailed steps or procedure of thisinvention:

IV. CONCLUSION AND ADVANTAGES

The invention proposes a data mining and least square estimationcombination framework to determine the overall systematic or biaserror(s) introduced by PMU and its instrumentation channel. The leastsquare estimator is built based on general TL PI model. In thisestimation, TL parameters are equally critical to the result assynchrophasor measurements. For most power systems, detailed TLspecifications can be found from their affiliated utility's EnergyManagement System (EMS). However, the accuracy level of EMS data is evenlower than PMU data itself, a sensitivity analysis we conductedshows: 1) TL parameter accuracy is essentially and linearly influentialto the state estimation result; 2) the more accurate the TL parameter isthe more disperse the erroneous synchrophasor measurements are from thegood synchrophasor measurements. These two features recall a clusteringalgorithm called density-based spatial clustering of applications withnoise (DBSCAN), this algorithm is capable of effectively filtering outthe true synchrophasor bias errors from the massive estimation resultdespite the inaccurate TL parameters. This invention presents a novelapproach for online calibration of PMU by using density-based spatialclustering. As compared to existing methods, this invention has twomajor advantages: 1) it identifies the overall bias errors introduced byboth PMU and its instrumentation channel; 2) it does not requireaccurate system model/parameters. Therefore, it is applicable across awide spectrum of practical conditions. In addition, one by-product ofthe proposed approach is more accurate TL impedance estimates forimproved system modeling, more accurate protective relay settings, andother related applications. Future more, this invention could: 1) extendthe invented framework to system level to achieve simultaneouscalibration of multiple PMUs; 2) decompose the spatial clusteringprocess so that state-of-the-art parallel computing techniques can beemployed to speed up the computation.

This invention has the following advantages:

-   -   1) Easier to implement, and cheaper in cost: No additional        measurement equipment is needed; Labors are saved since this        invention leads to an automatic calibration.    -   2) Faster operation: Online synchrophasor measurements data can        be directly imported into this calibration program and        monitoring the PMU accuracy.    -   3) Additional commercial value: The byproduct of this invention        is more accurate impedance parameters of the TL for EMS        database, which can benefit utilities on power flow modeling,        protective relay settings and other related work.

The framework of this bias error detection method is unique; itbasically combined the least square estimation and the clusteringalgorithm DBSCAN. The interactive of these two tools can effectivelydetect the bias errors.

The calibration process is also unique, it purely based on existing datalike PMU measurements and TL parameters from EMS, no additionalassumptions or calibration devices are needed. This method can solve twoproblems at one time. Beside bias error calibration, it also canidentify the true transmission line parameters.

V. APPENDIX I

Partial derivatives of impedance to the PMU measurements are presentedas follows:

$\begin{matrix}{\frac{\partial Z}{\partial V_{s}} = {{- \frac{2V_{s}e^{2i\; \theta_{Vs}^{\prime}}}{{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}}} - \frac{I_{r}{e^{i\; \theta_{Vs}^{\prime}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}} & (22) \\{\frac{\partial Z}{\partial V_{r}} = {\frac{2V_{r}e^{2i\; \theta_{Vr}^{\prime}}}{{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} + \frac{I_{s}{e^{i\; {({\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}} & (23) \\{\mspace{79mu} {\frac{\partial Z}{\partial I_{s}} = \frac{V_{r}{e^{i\; {({\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}} & (24) \\{\mspace{79mu} {\frac{\partial Z}{\partial I_{r}} = {- \frac{V_{s}{e^{i\; \theta_{Vs}^{\prime}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}}} & (25) \\{\frac{\partial Z}{{\partial\theta}\; V_{s}^{\prime}} = {{- \frac{2{iV}_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}{{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}}} - \frac{{iI}_{r}V_{s}{e^{i\; \theta_{Vs}^{\prime}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}} & (26) \\{\frac{\partial Z}{{\partial\theta}\; V_{r}^{\prime}} = {\frac{2{iV}_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}}{{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} - \frac{{iI}_{s}V_{r}{e^{i{({\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}} & (27) \\{\mspace{79mu} {{{\frac{\partial Z}{{\partial\theta}\; I_{s}^{\prime}} = \frac{{iI}_{s}V_{r}{e^{i\; {({\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}\left( {{V_{r}^{2}e^{2i\; \theta_{Vr}^{\prime}}} - {V_{s}^{2}e^{2i\; \theta_{Vs}^{\prime}}}} \right)}}{\left( {{I_{r}V_{s}e^{i\; \theta_{Vs}^{\prime}}} - {I_{s}V_{r}e^{i(\; {\theta_{Vr}^{\prime} + \theta_{Is}^{\prime}})}}} \right)^{2}}}\mspace{20mu} {A_{R} = {{real}\left( \frac{\partial Z}{\partial V_{s}} \right)}}},{A_{X} = {{imag}\left( \frac{\partial Z}{\partial V_{s}} \right)}},{B_{R} = {{real}\left( \frac{\partial R}{\partial V_{r}} \right)}},\mspace{20mu} {B_{X} = {{imag}\left( \frac{\partial R}{\partial V_{r}} \right)}},{C_{R} = {{real}\left( \frac{\partial Z}{\partial I_{s}} \right)}},{C_{X} = {{imag}\left( \frac{\partial Z}{\partial I_{s}} \right)}},\mspace{20mu} {D_{R} = {{real}\left( \frac{\partial Z}{\partial I_{r}} \right)}},{D_{X} = {{imag}\left( \frac{\partial Z}{\partial I_{r}} \right)}},{E_{R} = {{real}\left( \frac{\partial Z}{\partial\theta_{Vs}^{\prime}} \right)}},\mspace{20mu} {E_{X} = {{imag}\left( \frac{\partial Z}{\partial\theta_{Vs}^{\prime}} \right)}},{F_{R} = {{real}\left( \frac{\partial Z}{\partial\theta_{Vr}^{\prime}} \right)}},{F_{X} = {{imag}\left( \frac{\partial Z}{\partial\theta_{Vr}^{\prime}} \right)}},\mspace{20mu} {G_{R} = {{real}\left( \frac{\partial Z}{{\partial\theta}\; I_{s}^{\prime}} \right)}},{G_{X} = {{{imag}\left( \frac{\partial Z}{{\partial\theta}\; I_{s}^{\prime}} \right)}.}}}} & (28)\end{matrix}$

By the same means, the partial derivative equations of Y to each PMUcomponents can be generated and A_(B), B_(B) . . . G_(B) can becalculated accordingly. Due to the space limitation, they are notpresented here.

VI. APPENDIX II

A transmission line with two PMUs installed at both terminals issimulated in this study using Matlab/Simulink with specifications shownin Table VIII.

TABLE VIII SPECIFICATIONS OF THE SIMULATED TRANSMISSION LINE VariablesDescription, Unit Value R_(line) line resistance, Ω/km 0.013333 L_(line)line inductances, H/km 7.4342e−4 C_(line) line capacitances, F/km1.0001e−8 D length of line, km 150 f_(source) source frequency, Hz 60Voltage level kV 500

VII. REFERENCES

-   [1] D. Shi, D. J. Tylaysky, and N. Logic, “An Adaptive Method for    Detection and Correction of Errors in PMU Measurements,” IEEE Trans.    Smart Grid, vol. 3, no. 4, December 2012.-   [2] North American SynchroPhasor Initiative Performance & Standards    Task Team, “Synchrophasor Measurement Accuracy Characterization,”    November 2007.-   [3] North American SynchronPhasor Initiatibe (NASPI), Performance    and Standards Task Team (PSTT), “PMU System Testing and Calibration    Guide”, November 5, 2007.-   [4] IEEE Standard for Synchrophasor Measurements for Power Systems,    IEEE Std C37. 118. 1™-2011.-   [5] P. Komarnicki, C. Dzienis, Z. A. Styczynski, J. Blumschein,    and V. Centeno, “Practical Experience with PMU System Testing and    Calibration Requirements,” IEEE PES General Meeting, 2008.-   [6] K. Narendra, Z. Zhang, J. Lane, B. Lackey and E. Khan,    “Calibration and Testing of TESLA Phasor Measurement Unit (PMU)    Using Double F6150 Test Instrument,” Bulk Power System Dynamics and    Control—VII. Revitalizing Operational Reliability, 2007 iREP    Symposium, pp. 1-13, August 2007.-   [7] G. Stenbakken and T. Nelson, “Static Calibration and Dynamic    Characterization of PMUs at NIST,” IEEE PES General Meeting, pp.    1-4, June 2007.-   [8] G. Stenbakken and M. Zhou, “Dynamic Phasor Measurement Unit Test    System,” IEEE PES General Meeting, pp. 1-8, June 2007.-   [9] G. Stenbakken, T. Nelson, M. Zhou and V. Centeno, “Reference    Values for Dynamic Calibration of PMUs,” Proceedings of the 41st    Hawaii International Conference on System Sciences, pp. 1-6, January    2008.-   [10] Y. Hu and D. Novosel, “Progresses in PMU Testing and    Calibration,” in Proc. 3rd Int. Con. on Electric Utility    Deregulation and Restructuring and Power Technologies, Nanjing, pp.    150 155, 2008.-   [11] R. Pirret, “Testing and Calibration of Phasor Measurement    Units,” Technical note, Fluke Calibration, 2013.-   [12] L. Vanfretti, J. H. Chow, S. Sarawgi, and B. Fardanesh, “A    Phasor-data-based State Estimator Incorporating Phase Bias    Correction,” IEEE Trans. Power Syst., vol. 26, pp. 111-119, February    2011.-   [13] A. P. S. Meliopoulos, G. J. Cokkinides, F. Galvan, and B.    Fardanesh, “GPS-Synchronized Data Acquisition: Technology Assessment    and Research Issues,” Proceedings of the 39th Annual Hawaii    International Conference on System Science, Hawaii, January 4-7,    2006.-   [14] Q. Zhang, V. Vittal, G. Heydt, etc., “The Integrated    Calibration of Synchronized Phasor Measurement Data in Power    Transmission Systems,” IEEE Trans. Power Devliery, vol. 26, no. 4,    pp.2573-2581, 2011.-   [15] Z. Wu, L. Zora, and A. G. Phadke, “Simultaneous Transmission    Line Parameter and PMU Measurement Calibration,” IEEE PES General    Meeting, 2015.-   [16] D. Ritzmann, P. S. Wright, W. Holderbaum, etc., “A Method for    Accurate Transmission Line Impedance Parameter Estimation,” IEEE    Trans. Instru. and Meas., vol. 65, no. 10, October 2016.-   [17] P. Castello, M. Lixia, C. Muscas, and P. A. Pegoraro., “Impact    of the Model on the Accuracy of Synchrophasor Measurement,” IEEE    Trans. Instr. and Meass, vol. 61, no. 8, pp. 2179-2188, August 2012.-   [18] D. Shi, D. J. Tylaysky, K. Koellner, N. Logic, and D. E.    Wheeler, “Transmission Line Parameters Identification Using PMU    Measurements,” Eur. Trans. Elect. Power, vol. 21, no. 4, pp.    1574-1588, 2011.-   [19] C. S. Indulkar and K. Ramalingam, “Estimation of Transmission    Line Parameters from Measurements,” Int. J. Elect. Power Energy    Syst., vol. 30, no. 5, pp. 337-342, June 2008.-   [20] C. Borda, A. Olarte, and. H. Diaz, “PMU-based Line and    Transformer Parameter Estimation,” in Proc. IEEE/PES Power :Syst.    Conf Expo., March 2009, pp. 1-8.-   [21] K. Dasgupta and S. A. Soman, “Line Parameter Estimation Using    Phasor Measurements by the Total Least Squares Approach,” in Proc.    IEEE Power Energy Soc. General Meeting, 1-5, July 2013.-   [22] L. Philippot and J. C. Maun, “An Application of Synchronous    Phasor Measurement to the Estimation of the Parameters of an    Overhead Transmission Line,” in Proc. Conf. Fault Disturbance Anal.    Precise Meas. Power Syst., Arlington, Va., pp. 1-5, November 1995.-   [23] A. M. Dán and D. Raisz, “Estimation of Transmission Line    Parameters Using Wide-area Measurement Method,” Proc. IEEE Trondheim    PowerTech, pp. 1-6. June 2011.-   [24] J. Zaborszky and J. W. Rittenhouse, Electric Power    Transmission: The Power System in the Steady State. New York, N.Y.,    USA: Ronald, 1954.-   [25] Cauchy-Riernann Equations. [Online]. Available at:    http://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann equations.-   [26] D. Shi, Utilizing Synchrophasor Technology to Determine    Transmission Line Impedance Parameters, Arizona State University,    2009.-   [27] H. Zhou, X. Zhao, D. Shi, etc., “Calculating Sequence    Impedances of Transmission Line Using PMU Measurements,” 2015 IEEE    Power & Energy Society General Meeting, 2015.-   [28] M. Ester, H. P. Kriegel, J. Sander, and X. Xu, “A Density-Based    Algorithm for Discovering Clusters in Large Spatial Databases with    Noise,” Kdd., vol. 96, no. 34, 1996.

In summary, the present invention provides a novel approach for onlinecalibration of PMU by using density-based spatial clustering. Althoughthe potential methods and applications of using the same according tothe present invention have been described in the foregoing specificationwith considerable details, it is to be understood that modifications maybe made to the invention which do not exceed the claimed scope of theinvention and modified forms of the present invention done by othersskilled in the art to which the invention pertains will be consideredinfringements of this invention when those modified forms fall withinthe claimed scope of this invention.

What is claimed is:
 1. A method of online calibration of PhasorMeasurement Unit (PMU ) by using Density-Based Spatial Clustering ofApplications with Noise (DBSCAN) approach with relaxed assumptions, themethod comprising the following steps: 1) Collect the voltage phasors, V_(s) and V _(r), and current phasors, Ī_(s) and Ī_(r), from bothends/terminals (the sending end, denoted as s, and the receiving end,denoted as r) of a transmission line; 2) Calculate the transmission lineseries resistance R, series reactance X, and shunt susceptance B_(c), bysolving a set of nodal equations (3)-(4), whose solution is given byequations (9)-(11); $\begin{matrix}{{{\overset{\_}{I}}_{s} - {{\overset{\_}{V}}_{s} \cdot \frac{Y}{2}} + {\overset{\_}{I}}_{r} - {{\overset{\_}{V}}_{r} \cdot \frac{Y}{2}}} = 0} & (3) \\{{{\overset{\_}{V}}_{s} - {Z \cdot \left( {{\overset{\_}{I}}_{s} - {{\overset{\_}{V}}_{s} \cdot \frac{Y}{2}}} \right)} - {\overset{\_}{V}}_{r}} = 0} & (4) \\{R = {{real}\left( \frac{{V_{s}^{2}e^{i\; 2\; \theta_{V_{s}}^{\prime}}} - {V_{r}^{2}e^{i\; 2\; \theta_{V_{r}}^{\prime}}}}{{{I_{s} \cdot V_{r}}e^{i({\theta_{I_{s}}^{\prime} + \theta_{V_{r}}^{\prime}})}} - {{I_{r} \cdot V_{s}}e^{i\; \theta_{V_{s}}^{\prime}}}} \right)}} & (9) \\{X = {{imag}\left( \frac{{V_{s}^{2}e^{i\; 2\; \theta_{V_{s}}^{\prime}}} - {V_{r}^{2}e^{i\; 2\; \theta_{V_{r}}^{\prime}}}}{{{I_{s} \cdot V_{r}}e^{i({\theta_{I_{s}}^{\prime} + \theta_{V_{r}}^{\prime}})}} - {{I_{r} \cdot V_{s}}e^{i\; \theta_{V_{s}}^{\prime}}}} \right)}} & (10) \\{B_{c} = {{imag}\left( {2 \cdot \frac{{I_{s}e^{i\; \theta_{I_{s}}^{\prime}}} + I_{r}}{{V_{s}e^{i\; \theta_{V_{s}}^{\prime}}} + {V_{r}e^{i\; \theta_{V_{r}}^{\prime}}}}} \right)}} & (11)\end{matrix}$ 3) Obtain the reference values of series resistanceR_(EMS), series reactance X_(EMS), and shunt susceptance B_(EMS) fromthe database of the Energy Management System (EMS) and compare thecalculated parameters obtained from step 2) to these reference values toevaluate the difference between the calculated values (R, X, B_(c)) andthe referenced ones (R_(EMS), X_(EMS), B_(EMS)), which are denoted as∂R, ∂X, and ∂B_(c), where ∂R=R−R_(EMS), ∂X=X−X_(EMS),∂B_(c)=B_(c)=B_(EMS). 4) Use Equations (9)-(11) (which are derived basedon nodal equations for a transmission line) to be linearized in order toidentify the relationship between bias errors in the PMU measurements(∂V_(s), ∂V_(r), ∂I_(s), ∂I_(r), ∂θ′_(Vs), ∂θ′_(Vr), and ∂θ′_(Is)) andthe differences/errors in the calculated parameters, ∂R, ∂X, and ∂B_(c),as described by equations (12)-(14), or equation (15) in a matrixformat. With multiple sets of PMU measurements, equations (15) can beextended and finally solved with solution given by equation (19). Thatis, given one set of ∂R, ∂X, and ∂B_(c), the bias error in PMUmeasurements can be calculated accordingly.∂R=A _(R) ·∂V _(s) +B _(R) ·∂V _(r) +C _(R) ·∂I _(s) +D _(R) ·∂I _(r) +E_(R)··∂θ′_(V) _(s) +F _(R)·∂θ′_(V) _(r) +G _(R)·∂θ′_(I) _(s)   (12)∂X=A _(X) ·∂V _(s) +B _(X)·∂V _(r) +C _(X)·∂I_(s) +D _(X) ·∂I _(r) +E_(X)·∂θ′_(V) _(s) +F _(X)·∂θ′_(V) _(r) +G _(X)·∂θ′_(I) _(s)   (13)∂B _(c) =A _(B) ·∂V _(s) +B _(B) ·∂V _(r) +C _(B) ·∂I _(s) +D _(B) ·∂I_(r) +E _(B)·∂θ′_(V) _(s) +F _(X)·∂θ′_(V) _(r) +G _(X)·∂θ′_(I) _(s)  (14) $\begin{matrix}{\begin{bmatrix}{\partial R} \\{\partial X} \\{\partial B_{c}}\end{bmatrix} = {\begin{bmatrix}A_{R} & B_{R} & C_{R} & D_{R} & E_{R} & F_{R} & G_{R} \\A_{X} & B_{X} & C_{X} & D_{X} & E_{X} & F_{X} & G_{X} \\A_{B} & B_{B} & C_{B} & D_{B} & E_{B} & F_{B} & G_{B}\end{bmatrix} \cdot \begin{bmatrix}{\partial V_{s}} \\{\partial V_{r}} \\{\partial I_{s}} \\{\partial I_{r}} \\{\partial\theta_{V_{s}}^{\prime}} \\{\partial\theta_{V_{r}}^{\prime}} \\{\partial\theta_{I_{s}}^{\prime}}\end{bmatrix}}} & (15)\end{matrix}$F=(H ^(T) H)⁻¹ H ^(T) E  (19) 5) Add perturbance to reference values ofR, X, Bc, based on the reference values in the EMS database sincepractically the accurate reference values of R, X, and B_(c) can neverbe known. For example, if the reference value of R in the EMS databaseis R0, it is assumed that the actual R of the line lies within a certainrange of this R0, such as 0.8*R0˜1.2*R0, which is named as a feasibilityregion for the reference value of this parameter. 6) Create afeasibility region for each of these three referenced parametersR_(EMS), X_(EMS), and B_(EMS), and vary these referenced values withintheir corresponding feasibility region. Base on each set of newreference values to calculate ∂R, ∂X, and ∂B_(c) and calculate thecorresponding bias errors in the PMU measurements (∂V_(s), ∂V_(r),∂I_(s), ∂I_(r), ∂θ′_(Vs), ∂θ′_(Vr), and ∂θ′_(Is)) by solving equation(21), which is an extended version of equation (15).$\mspace{745mu} {{(21)\begin{bmatrix}{\partial V_{s}^{1}} & {\partial V_{s}^{2}} & \; & {\partial V_{s}^{M}} \\{\partial V_{r}^{1}} & {\partial V_{r}^{2}} & \; & {\partial V_{r}^{M}} \\{\partial I_{s}^{1}} & {\partial I_{s}^{2}} & \ldots & {\partial I_{s}^{M}} \\{\partial I_{r}^{1}} & {\partial I_{r}^{2}} & \ldots & {\partial I_{r}^{M}} \\{{\partial\theta}\; V_{s}^{\prime 1}} & {{\partial\theta}\; V_{s}^{\prime 2}} & \ldots & {{\partial\theta}\; V_{s}^{\prime \; M}} \\{{\partial\theta}\; V_{r}^{\prime 1}} & {{\partial\theta}\; V_{r}^{\prime 2}} & \; & {{\partial\theta}\; V_{r}^{\prime \; M}} \\{{\partial\theta}\; I_{s}^{\prime 1}} & {{\partial\theta}\; I_{s}^{\prime 2}} & \; & {{\partial\theta}\; I_{s}^{\prime \; M}}\end{bmatrix}} = {\left( {H^{T}H} \right)^{- 1}{H^{T}\begin{bmatrix}{\partial R_{1}^{1}} & {\partial R_{1}^{2}} & \; & {\partial R_{1}^{M}} \\{\partial X_{1}^{1}} & {\partial X_{1}^{2}} & \; & {\partial X_{1}^{M}} \\{\partial B_{c\; 1}^{1}} & {\partial B_{c\; 1}^{2}} & \; & {\partial B_{c\; 1}^{M}} \\{\partial R_{2}^{1}} & {\partial R_{2}^{2}} & \; & {\partial R_{2}^{M}} \\{\partial X_{2}^{1}} & {\partial X_{2}^{2}} & \ldots & {\partial X_{2}^{M}} \\{\partial B_{c\; 2}^{1}} & {\partial B_{c\; 2}^{2}} & \ldots & {\partial B_{c\; 2}^{M}} \\\vdots & \vdots & \ldots & \vdots \\\vdots & \vdots & \; & \vdots \\\vdots & \vdots & \; & \vdots \\{\partial R_{N}^{1}} & {\partial R_{N}^{2}} & \; & {\partial R_{N}^{M}} \\{\partial X_{N}^{1}} & {\partial X_{N}^{2}} & \; & {\partial X_{N}^{M}} \\{\partial B_{cN}^{1}} & {\partial B_{cN}^{2}} & \; & {\partial B_{cN}^{M}}\end{bmatrix}}}}$ 7) For each set of bias errors (∂V_(s), ∂V_(r),∂I_(s), ∂I_(r), ∂θ′_(Vs), ∂θ′_(Vr), and ∂θ′_(Is)) calculated in step 6),apply DB SCAN algorithm to cluster them. Each time using DBSCAN tocluster the set of bias errors, record the biggest cluster and theminimum searching distance. 8) Evaluate the clustering results for eachset of bias errors to find out a group of clusters C=[C₁, C₂, . . . ,C_(n)] with the biggest cluster size(s). 9) Check the searching distanceof all clusters of C to find out the cluster C_(m) with the minimumsearching distance. 10) The bias errors (∂V_(s), ∂V_(r), ∂I_(s), ∂I_(r),∂θ′_(Vs), ∂θ′_(Vr), and ∂θ′_(Is))to cluster C_(m) are the true biaserrors of the PMU measurements and the corresponding values of R_(EMS),X_(EMS), and B_(EMS) are the true value of the transmission lineparameters.
 2. The method of online calibration of Phasor MeasurementUnit (PMU) in claim 1, wherein the method identifies the overall biaserrors introduced by both PMU and its instrumentation channel, whereinthe method does not require accurate system model/parameters.
 3. Themethod of online calibration of Phasor Measurement Unit (PMU) in claim1, wherein the method is more accurate in transmission line (TL)impedance estimation for improved system modeling, and is more accuratein protective relay settings or other related applications.